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ABSTRACT 

In  tomographic  medical  devices  such  as  SPECT  or  PET 
cameras,  image  reconstruction  is  an  unstable  inverse  prob¬ 
lem,  due  to  the  presence  of  additive  noise.  A  new  family 
of  regularization  methods  for  reconstruction,  based  on  a 
thresholding  procedure  in  wavelet  and  wavelet  packet  de¬ 
compositions,  is  studied.  This  approach  is  based  on  the  fact 
that  the  decompositions  provide  a  near-diagonalization  of 
the  inverse  Radon  transform  and  of  the  prior  information 
on  medical  images.  Corresponding  algorithms  have  been 
developed  for  both  2-D  and  full  3-D  reconstruction.  These 
procedures  are  fast,  non-iterative,  flexible,  and  their  perfor¬ 
mances  outperform  Filtered  Back-Projection  and  iterative 
procedures  such  as  OS-EM. 


1.  INTRODUCTION 

We  are  interested  in  the  problem  of  tomographic  recon¬ 
struction  of  images  from  transmission  data,  which  we  call 
tomographic  projections  or  sinograms.  Although  the  work 
presented  here  has  a  wide  range  of  applications  for  vari¬ 
ous  tomographic  devices,  we  will  focus  on  medical  images 
with  SPECT  and  PET  cameras. 

A  section  of  the  object  observed  by  the  tomographic 
device  is  represented  by  a  2-D  discrete  image  /[ni,ri2]. 
An  estimation  of  /  must  be  computed  with  a  tomographic 
reconstruction  procedure  from  the  sinograms  produced  by 
tomographic  devices,  denoted  Y[t,  9 ],  and  defined  as: 

Y[t,e\  =  n(f[n1,n2])  +  W[t,e\  (i) 

where  {/[m, n2]}0<n1<jv1-i, o<n2<N2-i  is  the  observed 
image,  W  is  an  additive  noise,  and  1Z  is  the  discrete  Radon 
transform  which  models  the  tomographic  projection  pro¬ 
cess.  The  discrete  Radon  transform  is  derived  from  its  con¬ 
tinuous  version  1ZC,  which  is  equivalent  to  the  X-ray  trans¬ 
form  in  two  dimensions  and  is  defined  as  [1] 


(ftc/cXM) 


fc(xi,X2)S(xi  cos  8+x 2  sin  9—t)dx\dx 

(2) 


where  fc(x  1,2:2)  €E  L2(M2),  5  is  the  Dirac  mass,  6  £ 

[0,  27t)  ,  and  t  £  R.  There  are  several  different  ways  to 
define  the  discrete  Radon  transform  based  on  the  continu¬ 
ous  Radon  transform  [2],  Typically,  a  line  integral  along 
X\  cos  9  +  x 2  sin  9  =  t  is  approximated  by  a  summation 
of  the  pixel  values  inside  the  strip  t  —  1/2  <  ni  cos  9  + 
ri2  sin#  <  t+  1/2. 

The  noise  W  is  usually  modelled  as  a  Poisson,  or  some¬ 
times  Gaussian  noise.  However,  since  the  tomographic 
projections  Y  are  often  processed  to  incorporate  various 
corrections,  such  as  attenuation  correction,  scatter  correc¬ 
tion,  resolution  correction  or  geometric  correction,  the  re¬ 
sulting  noise  is  also  distorted  and  does  not  always  com¬ 
ply  with  such  prior  models.  The  present  approach  can  be 
adapted  to  different  types  of  noise,  including  the  case  when 
there  is  no  available  statistical  model  for  the  noise. 

A  tomographic  reconstruction  procedure  must  incorpo¬ 
rate  the  following  steps:  a  backprojection  may  be  viewed 
as  the  application  of  a  discretized  inverse  Radon  transform 
7Z~  1  on  the  tomographic  projections  Y .  This  can  be  di¬ 
rectly  computed  with  a  radial  interpolation  and  a  decon¬ 
volution  to  amplify  the  high  frequency  components  of  the 
tomographic  projections  Y  in  the  direction  of  t.  This  de- 
convolution  comes  from  the  fact  that  the  Radon  transform 
is  a  smoothing  transform.  Consequently,  backprojecting 
in  the  presence  of  additive  noise  is  an  ill-posed  inverse 
problem:  numerically  speaking,  a  direct  computation  of 
rRT1Y  =  f  +  Z  is  contaminated  by  a  large  additive  noise 
Z  =  7Z~1W,  which  means  that  a  regularization  has  to  be 
incorporated  in  the  reconstruction  procedure. 

Current  approaches  for  regularization  in  tomographic 
reconstruction  can  be  summarized  as  follows: 

1.  Filtered  Back-Projection  (FBP)  are  linear  filtering 
techniques  in  the  Fourier  space.  FBP  suffers  from 
performance  limitations  due  to  the  fact  that  the  si¬ 
nusoids  of  the  Fourier  basis  are  not  adapted  to  rep¬ 
resent  spatially  inhomogeneous  data  as  found  in  med¬ 
ical  images. 

2.  Iterative  statistical  model-based  techniques  are  de¬ 
signed  to  implement  Expectation-Maximization  (EM) 
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and  Maximum  A  Posteriori  (MAP)  estimators  [3, 

4].  In  some  cases,  these  approaches  can  provide  an 
improvement  over  FBP,  but  these  estimators  suffer 
from  the  following  drawbacks: 

•  Computation  time.  Almost  all  the  correspond¬ 
ing  algorithms  are  too  computer-intensive  for 
clinical  applications,  with  the  exception  of  OS- 
EM  [5]  (an  accelerated  implementation  of  an 
EM  estimator).  In  MAP  methods,  useful  pri¬ 
ors  usually  give  local  maxima,  and  the  com¬ 
putational  cost  of  relaxation  methods  is  pro¬ 
hibitive. 

•  Theoretical  understanding  and  justification. 

EM  estimators  lack  theoretical  foundations  to 
understand  and  characterize  the  estimation  er¬ 
ror.  Some  MAP  estimators  are  in  some  cases 
better  understood,  yet  no  optimality  for  a  re¬ 
alistic  model  has  been  established. 

•  Convergence.  EM  estimators  are  ill-conditioned, 
in  the  sense  that  the  corresponding  iterative 
algorithms  have  to  be  stopped  after  a  limited 
number  of  iterations.  Beyond  this  critical  num¬ 
ber,  the  noise  is  magnified,  and  EM  and  OS- 
EM  converge  to  a  non-ML  (Maximum  Like¬ 
lihood)  solution. 

This  study  aims  to  adress  these  limitations  by  building  a 
family  of  estimation  procedures  which  provide  better  nu¬ 
merical  results,  both  metrically  and  perceptually,  with  a 
sound  theoretical  basis  and  for  which  the  estimation  error  is 
understood  and  characterized.  The  estimators  should  also 
be  implemented  with  fast  and  flexible  algorithms. 

2.  THRESHOLDING  ESTIMATORS 

The  estimation  problem  in  (1)  is  also  equivalent  to  the  de- 
noising  problem 

X  =  f  +  Z  (3) 

where  X  =  R.  1  Y  and  Z  =  1Z~1W.  If  the  noise  Z  was 
white,  Donoho  and  Johnstone  have  established  [6]  that  a 
thresholding  estimator  in  a  properly  selected  vector  family 

®  =  {5,ni,7i2)5,n1,n2}o<ni<JVi-i,o<n2<JV2-i.  typically  a 
wavelet  basis,  would  be  optimal  to  recover  spatially  inho¬ 
mogeneous  data  as  found  in  medical  images.  A  threshold¬ 
ing  estimator  F  of  f  in  B  is  defined  as 

^  ^  '  Pri  \  .r/,2  ((X,  @711,122)}  9ni,ri2  ’ 

111  ,12 2 

where  pni,n2  is  a  thresholding  operator.  In  our  situation, 
the  choice  of  B  does  not  only  depend  on  the  prior  informa¬ 
tion  on  the  object  /,  but  also  on  the  noise  Z,  whose  behav¬ 
ior  is  very  specific  due  to  the  fact  that  a  backprojection  has 
been  applied  on  it. 


The  assumption  underlying  thresholding  estimators  is 
that  each  coefficient  in  the  decomposition  B  can  be  esti¬ 
mated  independently  without  a  loss  of  performance.  As 
a  consequence,  such  estimators  are  efficient  if  the  coeffi¬ 
cients  of  the  noise  and  of  the  object  to  be  recovered  are 
indeed  nearly  independent  in  B.  This  means  that  B  must 
provide  a  near-diagonalization  of  the  noise  Z  and  of  the 
prior  information  on  the  image  /. 

The  image  /  is  a  spatially  inhomogeneous,  piece-wise 
regular  signal,  which  is  compactly  represented  in  a  wavelet 
decomposition.  To  obtain  a  diagonal  representation  of  the 
noise  Z,  we  want  to  find  a  decomposition  in  which  the  in¬ 
verse  Radon  transform  is  nearly  diagonal.  Since  the  inverse 
Radon  transform  is  a  Calderon-Zygmund  operator  [7],  it  is 
also  nearly-diagonal  in  a  wavelet  basis. 

These  two  properties  of  wavelet  bases  led  Donoho  [8] 
to  suggest  the  use  of  thresholding  estimators  in  wavelet 
bases  for  several  linear  inverse  problems,  including  the  in¬ 
version  of  the  Radon  transform.  Donoho  established  the 
aymptotic  optimality  (minimax  sense)  of  this  approach,  called 
Wavelet- Vaguelette  Decomposition  (WVD),  and  its  supe¬ 
riority  with  respect  to  other  approaches  such  as  Filtered 
BackProjection.  However,  the  WVD  as  studied  by  Donoho 
is  a  theoretical  concept  developped  for  continuous  mod¬ 
els,  based  on  the  assumption  that  the  additive  noise  W  is 
Gaussian  white,  and,  despite  numerical  implementations 
and  refinements  by  other  researchers [9,  10],  this  technique 
is  not  as  such  adapted  to  the  tomographic  reconstruction  of 
real  medical  images.  In  the  meanwhile,  Kalifa  and  Mallat 
[11]  generalized  Donoho’s  approach  to  adapt  it  to  different 
types  of  decompositions,  not  restricted  to  wavelet  bases. 

Wavelet  packet  bases  are  other  decompositions  which 
can  provide  a  compact  representation  of  the  observed  im¬ 
age  /,  and  since  a  wavelet  packet  transform  provides  a 
more  accurate  segmentation  of  the  frequency  domain  than 
a  wavelet  transform,  this  improves  the  near-diagonalization 
of  the  noise  Z.  Besides,  a  wavelet  packet  basis  can  be  adap¬ 
tively  chosen  from  a  dictionary  of  different  wavelet  packet 
bases.  We  shall  see  in  the  next  section  that  this  enables  us 
to  optimize  the  choice  of  the  wavelet  packet  transform  to 
the  specific  type  of  observed  image  and  to  the  nature  of  the 
backprojected  noise  Z. 

3.  CHOICE  OF  WAVELET  PACKET 
DECOMPOSITION 

Suppose  we  have  a  dictionary  {B': }7  of  orthogonal  wavelet 
packet  bases.  The  empirical  best  basis  B11  for  estimating  / 
is  obtained  by  minimizing  an  estimation  of  the  final  estima¬ 
tion  error  (risk),  using  the  best  basis  algorithm  of  Coifman 
and  Wickerhauser  [12],  For  example,  the  final  estimation 
error  can  be  estimated  with  a  Stein  Unbiased  Risk  Estima¬ 
tor  (SURE)  [13].  Another  possibility  is  to  take  advantage 
of  the  fact  that  phantom  images  are  usually  available  in  to- 


mography  to  model  the  observed  images,  and  can  be  used 
as  reference  images  to  minimize  the  ideal  projection  risk, 
as  defined  in  [6]. 

In  both  cases,  a  numerical  model  of  the  noise  Z  is  com¬ 
puted.  Depending  on  the  fact  that  the  original  noise  W  is 
modelled  as  Gaussian  or  Poisson  noise,  or  not  modelled  by 
a  statistical  prior,  different  procedures  can  be  used.  The  nu¬ 
merical  model  of  the  noise  Z  is  taken  as  input  when  com¬ 
puting  the  choice  of  the  best  basis,  which  is  specifically 
adapted  not  only  to  the  nature  of  the  data  but  also  to  the 
nature  of  the  noise  Z.  This  adaptive  approach  enables  us  to 
improve  greatly  the  performance  of  the  decomposition,  as 
compared  to  a  classical  wavelet  basis. 

4.  RECONSTRUCTION  ALGORITHM  AND 
NUMERICAL  RESULTS 

The  tomographic  reconstruction  algorithm  is  decomposed 
in  the  following  steps: 

1 .  Backprojection  without  regularization  of  the  tomo¬ 
graphic  projections  Y  to  obtain  the  backprojected 
image  X  =  f  +  Z. 

2.  (Optional)  Computation  of  the  best  wavelet  packet 
basis  B 71  optimized  for  the  specific  image  to  be  re¬ 
stored.  The  best  basis  can  be  recomputed  for  each 
image,  or  can  have  been  computed  once  and  stored 
in  advance,  to  save  computation  time.  However  the 
computation  time  of  the  best  basis  algorithm  is  very 
short. 

3.  Wavelet  packet  transform  of  the  backprojected  im¬ 
age  X  in  the  best  basis  B 71  to  obtain  the  wavelet 
packet  coefficients  {(X,  gni,n2  )}  ni,n2* 

4.  Thresholding  on  the  wavelet  packet  coefficients. 

5.  Inverse  wavelet  packet  transform  of  the  thresholded 
coefficients  to  obtain  the  estimated  image  F. 

The  wavelet  packet  transform  and  its  inverse  are  computed 
with  fast  filter  bank  algorithms  of  complexity  O(N)  for 
signals  of  N  samples[14]. 

Figure  1  compares  preliminary  numerical  results  com¬ 
puted  on  SPECT  clinical  data,  using  an  OS-EM  reconstruc¬ 
tion,  a  Filtered  Back-Projection  and  a  wavelet-packet  based 
reconstruction.  The  OS-EM  reconstructed  image  is  very 
smooth  because  the  OS-EM  algorithm  has  to  be  stopped 
after  a  limited  number  of  iterations,  otherwise  the  noise  is 
strongly  amplified  and  the  algorithm  converges  to  a  very 
noisy  reconstructed  image.  On  the  other  hand,  the  FBP- 
reconstructed  image  is  corrupted  by  a  significant  noise  and 
artifacts,  which  cannot  be  reduced  unless  the  reconstructed 
image  becomes  extremely  smoothed.  With  the  wavelet  packet 
reconstruction  algorithms,  the  amount  of  smoothness  of  a 
reconstructed  image  can  be  controlled  precisely,  while  the 
noise  is  reduced  significantly  as  compared  to  an  image  re¬ 
constructed  with  FBP  or  OS-EM.  Besides,  wavelet  packet 


reconstruction  methods  are  adapted  to  a  broader  range  of 
images  (different  objects,  low  or  high  count)  than  OS-EM 
and  FBP.  The  superiority  of  the  wavelet-packet  based  re¬ 
construction  method  over  FBP  and  OS-EM  is  currently  be¬ 
ing  established  through  ROC  studies  on  various  SPECT 
and  PET  data. 


5.  EXTENSION  TO  3-D  RECONSTRUCTION 

So  far,  the  wavelet  packet  reconstruction  has  been  presented 
for  2-D  reconstruction  of  slices.  We  now  consider  that  we 
have  a  series  of  tomographic  projections  of  N3  translated 
2-D  slices  of  the  observed  object,  i.e.,  that  we  have  3D 
data.  When  necessary,  the  tomographic  projections  have 
been  transformed  via  rebinning  techniques  in  order  to  ob¬ 
tain  tomographic  projections  of  2-D  slices:  this  approach  is 
not  necessary  for  SPECT  images,  but  is  increasingly  com¬ 
mon  in  3-D  PET  imaging.  It  is  useful  to  take  advantage 
of  the  correlations  of  the  signal  in  the  transaxial  direction 
to  obtain  a  better  discrimination  between  the  information 
and  the  noise.  In  this  case,  a  regularization  is  computed  on 
the  whole  3-D  data,  but  the  backprojections  are  still  com¬ 
puted  slice  by  slice,  since  the  Radon  transform  is  still  only 
a  bidimensional  operator. 

The  3-D  thresholding  reconstruction  algorithm  is  sum¬ 
marized  as  follows.  Details  will  be  provided  at  the  confer¬ 
ence. 

•  Each  2-D  slice  is  reconstructed  and  partly  regular¬ 
ized  using  the  wavelet  packet  reconstruction  algo¬ 
rithm.  The  regularization  is  mild  in  order  to  preserve 
as  much  information  as  possible. 

•  The  remaining  noise  (after  the  2-D  regularization)  is 
nearly  diagonalized  in  a  wavelet  decomposition.  To 
take  fully  advantage  of  the  3-D  information  in  the 
data,  a  3-D  dyadic  wavelet  transform  is  applied  on 
the  3-D  data,  with  wavelets  adaptively  oriented  per¬ 
pendicularly  to  the  singularities  of  the  signal.  This 
directional  selectivity  enables  us  to  maximize  the  cor¬ 
relation  between  the  vectors  of  the  wavelet  family 
and  the  information  of  the  signal.  The  efficiency  of 
noise  removal  is  thus  greatly  improved. 
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